## Chapter 5 (COVID Analysis) 
## Winter 2023  

## create an indicator for party*wave 
data <- data %>% 
  mutate(party_year = case_when(
    dem == 1 & year == 2 ~ 1, 
    dem == 1 & year == 3 ~ 2, 
    dem == 1 & year == 4 ~ 3, 
    dem == 0 & year == 2 ~ 4, 
    dem == 0 & year == 3 ~ 5, 
    dem == 0 & year == 4 ~ 6))

#####################################
## Figure 5.1: Behaviors Over Time ##
#####################################

p <- data %>% 
        filter(year > 1) %>% 
        ggplot(aes(apstable, behavior,
                   group = rep)) + 
        geom_smooth(aes(color = as.factor(rep),
                        fill = as.factor(rep), 
                        alpha = 0.65), se = T) + 
        scale_color_manual(name = "",
                           values = c("black","grey70"), 
                           labels = c("Dem","Rep")) + 
        scale_fill_manual(name = "",
                          values = c("black","grey70"), 
                          labels = c("Dem","Rep")) + 
        scale_x_continuous(breaks = ap_quantiles, 
                           labels = c("5th \n pct","10th \n pct","33rd \n pct","50th \n pct","67th \n pct", "90th \n pct","95th \n pct")) + 
        xlab("Animus") +
        ylab("Number of Behaviors (0-14)") + 
        #ggtitle("Preventative Behaviors") + 
        labs(color = "") + 
        facet_grid(cols = vars(factor_year)) + 
        theme_bw() + 
        theme(plot.title = element_text(hjust=0.5),
              #aspect.ratio = 1, 
              legend.position = "none", 
              panel.grid = element_blank())    

## Preventative Behaviors, all 3 waves 
pdf(file = "Figures/Druckman_etal_Figure_5_1.pdf", 
    height = 8.5,
    width = 11) 
print(p)
dev.off() 

## extract values at 95% percentile  
## this takes the highest x value for Ds minus Rs in each panel 
behavior_data <- ggplot_build(p)$data[[1]] 
behavior_data$y[80] - behavior_data$y[160]  ## 1.3 (April 2020, this is panel 1) 
behavior_data$y[240] - behavior_data$y[320]  ## 2.05 (October 2020, this is panel 2) 
behavior_data$y[400] - behavior_data$y[480]  ## 3.56 (April 2021, this is panel 3) 


##############################
## Figure 5.2: Mask Wearing ## 
##############################

pdf(file = "Figures/Druckman_etal_figure_5_2.pdf", 
    height = 8.5,
    width = 11) 
print(data %>% 
        filter(year > 1) %>% 
        ggplot(aes(apstable, facemask,
                   group = rep)) + 
        geom_smooth(aes(color = as.factor(rep),
                        fill = as.factor(rep)),se = T) + 
        scale_color_manual(name = "",
                           values = c("black","grey70"), 
                           labels = c("Dem","Rep")) + 
        scale_fill_manual(name = "",
                           values = c("black","grey70"), 
                           labels = c("Dem","Rep")) + 
        scale_x_continuous(breaks = ap_quantiles, 
                           labels = c("5th \n pct","10th \n pct","33rd \n pct","50th \n pct","67th \n pct", "90th \n pct","95th \n pct")) + 
        xlab("Animus") +
        ylab("Likelihood of Mask Wearing") + 
        ylim(c(0,1)) + 
        #ggtitle("Mask Wearing over Time") + 
        labs(color = "") + 
        facet_grid(cols = vars(factor_year)) + 
        theme_bw() + 
        theme(plot.title = element_text(hjust=0.5),
              legend.position = "none", 
              #legend.box = "vertical",
              panel.grid = element_blank())) 
dev.off() 

## fraction wearing a mask over time 
data %>% 
  group_by(year) %>% 
  summarise(msk = mean(facemask,na.rm=T))

## is there a significant increase in the October 2020 data? 
library(mgcv) 
summary(gam(facemask ~ s(apstable),
            subset = dem == 1 & year == 3,
            data = data)) 
## not quite 

#####################################################
## Figure 5.3: Likelihood of Receiving the Vaccine ##
##################################################### 

## Will you get vaccinated (wave 4) 
## rescale so that already vaccinated & extremely likely is the top category 
data <- data %>% 
  mutate(vax_prob = case_when(
    w7vaccine < 4 ~ 4,
    w7vaccine == 4 ~ 3, 
    w7vaccine == 5 ~ 2, 
    w7vaccine == 6 ~ 1))  
data$been_vaccinated <- ifelse(data$w7vaccine==1,1,0) 

pdf(file = "Figures/Druckman_etal_Figure_5_3.pdf", 
    height = 11,
    width = 8.5) 
print(data %>% 
  filter(year == 4) %>% 
  ggplot(aes(apstable, vax_prob,
             group = dem)) + 
  geom_smooth(aes(color = as.factor(dem),
                  fill = as.factor(dem)), se = T) + 
  scale_color_manual(values = c("grey70","black"), 
                     labels = c("Republicans","Democrats")) + 
  scale_fill_manual(values = c("grey70","black"), 
                       labels = c("Republicans","Democrats")) + 
  scale_x_continuous(breaks = ap_quantiles, 
                     labels = c("5th \n pct","10th \n pct","33rd \n pct","50th \n pct","67th \n pct", "90th \n pct","95th \n pct")) + 
  xlab("Animus") +
  ylab("Likelihood of Vaccination") + 
  scale_y_continuous(breaks = c(1:4),
                     limits = c(1,4),
                     labels = c("Not at All \n Likely","Not Too \n Likely",
                                "Somewhat \n Likely", "Extremely Likely \n / Already Vaccinated")) + 
  #ggtitle("Likelihood of Vaccination") + 
  labs(color = "") + 
  theme_bw() + 
  theme(plot.title = element_text(hjust=0.5),
        legend.position = "none", 
        panel.grid = element_blank())) 
dev.off() 

######################################
## Figure 5.4: Safety of Behaviors  ##
###################################### 

## faceted by year 
pdf(file = "Figures/Druckman_etal_Figure_5_4.pdf", 
    height = 8.5,
    width = 11) 
## Safety (Waves 3 & 4 only)
data$safe <- (-1*data$unsafe) + 5 
print(data %>% 
        filter(year > 2) %>% 
        ## wave 4: no vaccine (for comparison) 
        #filter(safeexpcond1novac2vac==1) %>% 
        ggplot(aes(apstable, safe,
                   group = rep)) + 
        geom_smooth(aes(color = as.factor(rep),
                        fill = as.factor(rep)), se = T) + 
        scale_color_manual(name = "",
                           values = c("black","grey70"), 
                           labels = c("Democrat","Republican")) +
        scale_fill_manual(name = "",
                           values = c("black","grey70"), 
                           labels = c("Democrat","Republican")) +
        scale_x_continuous(breaks = ap_quantiles, 
                           labels = c("5th \n pct","10th \n pct","33rd \n pct","50th \n pct","67th \n pct", "90th \n pct","95th \n pct")) + 
        scale_y_continuous(breaks = c(1:4),
                           limits = c(1,4),
                           labels = c("Not at All \n Safe","Not Too \n Safe",
                                      "Somewhat \n Safe", "Very \n Safe")) + 
        xlab("Animus") +
        ylab("Perceived Safety of Various Activities") + 
        #ylim(c(1,4)) + 
        #ggtitle("Perceived Safety of Various Activities") + 
        labs(color = "", fill="") +
        facet_grid(cols = vars(factor_year)) + 
        theme_bw() + 
        theme(plot.title = element_text(hjust=0.5),
              legend.position = "none", 
              panel.grid = element_blank()))  
dev.off() 


## regression model 
## make exp. treatment 0/1 so it's more interpretable in the regression model 
data$vac_treat <- ifelse(data$safeexpcond1novac2vac==2,1,
                         ifelse(data$safeexpcond1novac2vac==1,0,NA)) 
summary(lm(safe ~ apstable*rep*vac_treat,
           data = data, 
           subset = year == 4)) 
summary(lm(safe ~ apstable*rep*vac_treat*been_vaccinated,
           data = data, 
           subset = year == 4)) 
summary(lm(safe ~ rep*vac_treat,
           data = data, 
           subset = year == 4)) 
## no effects 
## even on avg, no effect of vaccinated treatment 

## just look at partisan means 
data %>% 
  filter(year==4 & !is.na(vac_treat)) %>% 
  group_by(vac_treat,rep) %>% 
  summarise(safetybar = mean(safety,na.rm=T)) 
## diff-in-diff is basically 0 

## do vaccinated people think activities are safer in W4 than they did in W3? 
data <- data %>% 
  group_by(respondent_id) %>% 
  mutate(safe_lag = dplyr::lag(safety, n = 1, order_by = year)) %>% 
  rowwise(.) %>% 
  mutate(delta_safety = safety - safe_lag) %>% 
  ungroup(.) 

summary(lm(delta_safety ~ been_vaccinated, 
           data = data))
summary(lm(delta_safety ~ been_vaccinated*apstable*rep, 
           data = data))
  

#################################
## Figure 5.5: Policy Support  ## 
#################################

pdf(file = "Figures/Druckman_etal_Figure_5_5.pdf",
    height = 8.5,
    width =  11) 
print(data %>% 
  filter(year >2) %>% 
  ggplot(aes(apstable, policy,
             group = rep)) + 
  geom_smooth(aes(color = as.factor(rep),
                  fill = as.factor(rep)), se = T) + 
  scale_color_manual(name = "",
                     values = c("black","grey70"), 
                     labels = c("Democrat","Republican")) +
  scale_fill_manual(name = "",
                     values = c("black","grey70"), 
                     labels = c("Democrat","Republican")) +
  scale_x_continuous(breaks = ap_quantiles, 
                     labels = c("5th \n pct","10th \n pct","33rd \n pct","50th \n pct","67th \n pct", "90th \n pct","95th \n pct")) + 
  scale_y_continuous(breaks = c(1:4),
                      limits = c(1,4),
                      labels = c("Strongly Disagree","Somewhat Disagree",
                                  "Somewhat Agree", "Strongly Agree")) +   
  xlab("Animus") +
  ylab("Policy Position") + 
  #ylim(c(0,1)) + 
  #ggtitle("Support for COVID-19 Related Policies") + 
  labs(color = "") + 
  facet_grid(cols = vars(factor_year)) + 
  theme_bw() + 
  theme(plot.title = element_text(hjust=0.5),
        legend.position = "none", 
        panel.grid = element_blank())) 
dev.off()


#####################################
## Figure 5.7: Moderation & Safety ##
#####################################

## group severity into two bins, too much jumpiness with all four categories 
data$cm <- ifelse(data$covidmax == 1 | data$covidmax == 2,1,
                  ifelse(data$covidmax==3 | data$covidmax == 4,2,NA)) 
## put everyone into one plot 
data <- data %>% 
  mutate(cm_party = case_when( 
    cm == 1 & dem == 0 ~ 1, 
    cm == 2 & dem == 0 ~ 2, 
    cm == 1 & dem == 1 ~ 3, 
    cm == 2 & dem == 1 ~ 4)) 
##  make a factor variable for the labels 
data$sev_label <- factor(data$cm, 
                         labels = c("Low Consequences",
                                    "High Consequences"))

## plot, pool waves 3 and 4 
pdf(file = "Figures/Druckman_etal_figure_5_7.pdf", 
    height = 8.5,
    width = 11) 
print(data %>% 
        filter(!is.na(cm) & year > 2) %>%  
        ggplot(aes(apstable,safety),
               group = dem) + 
        geom_smooth(aes(color = as.factor(dem),
                        fill = as.factor(dem)), 
                    method = "gam", se = T) + 
        scale_color_manual(name = "",
                           values = c("grey70","black"), 
                           labels = c("Republican","Democrat")) + 
        scale_fill_manual(name = "",
                           values = c("grey70","black"), 
                           labels = c("Republican","Democrat")) + 
        scale_x_continuous(breaks = ap_quantiles, 
                           labels = c("5th \n pct","10th \n pct","33rd \n pct","50th \n pct","67th \n pct", "90th \n pct","95th \n pct")) + 
        scale_y_continuous(breaks = c(1:4),
                           limits = c(1,4),
                           labels = c("Not at All \n Safe","Not Too \n Safe",
                                      "Somewhat \n Safe", "Very \n Safe")) + 
        #ylim(c(1,4))+
        xlab("Animus") +
        ylab("Perceived Safety of Various Activities") + 
        #ggtitle("Perceived Safety & COVID-19 Severity") + 
        labs(color = "") + 
        facet_grid(cols = vars(sev_label)) + 
        theme_bw() + 
        theme(plot.title = element_text(hjust=0.5),
              legend.position = "none", 
              panel.grid = element_blank()))  
dev.off()

##########################################
## Figure 5.6: Moderation and Behaviors ##
########################################## 

## pool waves 3 and 4 
pdf(file = "Figures/Druckman_etal_figure_5_6.pdf", 
    height = 8.5,
    width = 11) 
print(data %>% 
        filter(!is.na(sev_label) & year > 2) %>%  
        ggplot(aes(apstable,behavior),
               group = dem) + 
        geom_smooth(aes(color = as.factor(dem),
                        fill = as.factor(dem)), 
                    method = "gam", se = T) + 
        scale_color_manual(name = "",
                           values = c("grey70","black"), 
                           labels = c("Republican","Democrat")) + 
        scale_fill_manual(name = "",
                           values = c("grey70","black"), 
                           labels = c("Republican","Democrat")) + 
        scale_x_continuous(breaks = ap_quantiles, 
                           labels = c("5th \n pct","10th \n pct","33rd \n pct","50th \n pct","67th \n pct", "90th \n pct","95th \n pct")) + 
        ylim(c(1,9))+
        xlab("Animus") +
        ylab("Number of Behaviors (0-14)") + 
        #ggtitle("Preventative Behaviors & COVID-19 Severity") + 
        labs(color = "") + 
        facet_grid(cols = vars(sev_label)) + 
        theme_bw() + 
        theme(plot.title = element_text(hjust=0.5),
              legend.position = "none", 
              panel.grid =element_blank()))
dev.off()

#############################################
## Figure 5.8: Policy Support and Salience ##
#############################################

## Policy Scale 
pdf(file = "Figures/Druckman_etal_Figure_5_8.pdf", 
    height = 8.5,
    width = 11) 
print(data %>% 
        filter(!is.na(cm) & year > 2) %>%  
        ggplot(aes(apstable,policy),
               group = dem) + 
        geom_smooth(aes(color = as.factor(dem),
                        fill = as.factor(dem)), 
                    method = "gam", se = T) + 
        scale_color_manual(name = "",
                           values = c("grey70","black"), 
                           labels = c("Republican","Democrat")) + 
        scale_fill_manual(name = "",
                          values = c("grey70","black"), 
                          labels = c("Republican","Democrat")) + 
        scale_x_continuous(breaks = ap_quantiles, 
                           labels = c("5th \n pct","10th \n pct","33rd \n pct","50th \n pct","67th \n pct", "90th \n pct","95th \n pct")) + 
        #ylim(c(1,4))+
        scale_y_continuous(breaks = c(1:4),
                           limits = c(1,4),
                           labels = c("Strongly Disagree","Somewhat Disagree",
                                      "Somewhat Agree", "Strongly Agree")) +  
        xlab("Animus") +
        ylab("Policy Position") + 
        #ggtitle("Policy Support & COVID-19 Severity") + 
        labs(color = "") + 
        facet_grid(cols = vars(sev_label)) + 
        theme_bw() + 
        theme(plot.title = element_text(hjust=0.5),
              legend.position = "none", 
              panel.grid = element_blank()))  
dev.off()


#######################################
## Table 5.1: First Difference Model ##
####################################### 

## Look at people who get worse experiences over time: do they do more protect activities? 
## For this, you need models, not a graph 

## First, calculate the change in the DVs and the covidmax variable 
data <- data %>% 
  group_by(respondent_id) %>% 
  mutate(cm_lag = dplyr::lag(covidmax, n = 1, order_by = year), 
         safe_lag = dplyr::lag(safety, n = 1, order_by = year), 
         behavior_lag = dplyr::lag(behavior, n = 1, order_by = year),
         policy_lag = dplyr::lag(policy, n = 1 , order_by = year)) %>% 
  rowwise(.) %>% 
  mutate(delta_cm = covidmax - cm_lag, 
         delta_safety = safety - safe_lag, 
         delta_behavior = behavior - behavior_lag, 
         delta_policy = policy - policy_lag) %>% 
  ungroup(.) 

## Worse/Better Cases 
100*prop.table(table(data$delta_cm))
## 29% imagine a less severe case, 50% no change, 21% imagine worse 

## Run the models 
## Pooling both parties                                   
m1 <- lm(delta_safety ~ delta_cm, data = data) 
m2 <- lm(delta_behavior ~ delta_cm, data = data)
m3 <- lm(delta_policy ~ delta_cm, data = data)
## Just look among Reps 
m4 <- lm(delta_safety ~ delta_cm, subset = rep == 1, data = data) 
m5 <- lm(delta_behavior ~ delta_cm, subset = rep == 1, data = data)
m6 <- lm(delta_policy ~ delta_cm, subset = rep == 1, data = data)

## export to a table 
models <- list()
models[["Safety \n (All)"]] <- m1 
models[["Behavior \n (All)"]] <- m2 
models[["Policy \n (All)"]] <- m3 
models[["Safety \n (Rs Only)"]] <- m4 
models[["Behavior \n (Rs Only)"]] <- m5
models[["Policy \n (Rs Only)"]] <- m6 

## clean up the GOF statistics (just add N/R2)
gm <- gof_map
gm$omit <- T 
gm$omit[2] <- F
gm$omit[5] <- F 

msummary(models,
         coef_map = c("(Intercept)" = "Intercept",
                      "delta_cm" = "Change in COVID Severity"), 
         stars = T, 
         #gof_map = gm, 
         #escape = F, 
         output =  "Tables/raw/Table_5_1.html")

## As a panel model (robustness check)
summary(felm(safety ~ covidmax | year + uniqueid,
             data = data))
summary(felm(behavior ~ covidmax | year + uniqueid,
             data = data))
summary(felm(policy ~ covidmax | year + uniqueid,
             data = data))
## same results: perceive less safety, more behaviors, and more policy support 

###################################
## Panel Models for the Appendix ##
################################### 

## Here, all models w/ just the April 2020 data are in the NHB Paper 
## We just need the panel models here 

m1 <- felm(safety ~ year*apstable*rep | uniqueid + year, 
           subset = partyid_w1 != 4 & partychange == 0,
             data = data)
m2 <- felm(behavior ~ year*apstable*rep | uniqueid + year, 
           subset = partyid_w1 != 4 & partychange == 0,
             data = data)
## policies change in W3, so only look at waves 3-4 
m3 <- felm(policy ~ year*apstable*rep | uniqueid + year, 
           subset = partyid_w1 != 4 & partychange == 0 & year > 2,
             data = data)  
m4 <- felm(facemask ~ year*apstable*rep | uniqueid + year, 
           subset = partyid_w1 != 4 & partychange == 0,
           data = data)

## output to tables 
## output this to a table 
models <- list()
models[["Behavior"]] <- m2 
models[["Face Mask"]] <- m4 
models[["Safety"]] <- m1  
models[["Policy"]] <- m3 

## add the N 
rows <- tribble(~term, ~"m1",~"m2",~"m3",~"m4",
                'N',m1$N,m2$N,m3$N,m4$N)
attr(rows,'position') <- 15 

## clean up the GOF output (just N & R-Squared)
gm <- modelsummary::gof_map 
gm$omit <- TRUE 
gm$omit[5] <- FALSE 

msummary(models,
         add_rows = rows, 
         coef_map = c("apstable" = "Animosity",
                      "rep" = "Republican",
                      "year" = "Wave",
                      "year:rep" = "Republican*Wave",
                      "year:apstable" = "Animosity*Wave",
                      "apstable:rep" = "Animosity*Republican", 
                      "year:apstable:rep" = "Animosity*Republican*Wave"), 
         gof_map = gm, 
         stars = T, 
         output =  "Tables/raw/Table_A5_1.html")  

## Can we exaplin away these effects? 
## Run models controlling for partisan strength & partisan social identity 

## estimate models with partisan strength 
data$pid_strength <- abs(data$partyid_w1 - 4) 
m1a <- felm(safety ~ year*apstable*rep + year*pid_strength*rep | uniqueid + year, 
           subset = partyid_w1 != 4 & partychange == 0,
           data = data)
m2a <- felm(behavior ~ year*apstable*rep + year*pid_strength*rep | uniqueid + year, 
           subset = partyid_w1 != 4 & partychange == 0,
           data = data)
## policies change in W3, so only look at waves 3-4 
m3a <- felm(policy ~ year*apstable*rep + year*pid_strength*rep | uniqueid + year, 
           subset = partyid_w1 != 4 & partychange == 0 & year > 2,
           data = data)  
m4a <- felm(facemask ~ year*apstable*rep + year*pid_strength*rep | uniqueid + year, 
           subset = partyid_w1 != 4 & partychange == 0,
           data = data) 


## also add in partisan social identity 
m1b <- felm(safety ~ year*apstable*rep + year*pid_strength*rep + year*sociden_pid*rep | uniqueid + year, 
            subset = partyid_w1 != 4 & partychange == 0,
            data = data)
m2b <- felm(behavior ~ year*apstable*rep + year*pid_strength*rep + year*sociden_pid*rep | uniqueid + year, 
            subset = partyid_w1 != 4 & partychange == 0,
            data = data)
## policies change in W3, so only look at waves 3-4 
m3b <- felm(policy ~ year*apstable*rep + year*pid_strength*rep + year*sociden_pid*rep | uniqueid + year, 
            subset = partyid_w1 != 4 & partychange == 0 & year > 2,
            data = data)  
m4b <- felm(facemask ~ year*apstable*rep + year*pid_strength*rep + year*sociden_pid*rep | uniqueid + year, 
            subset = partyid_w1 != 4 & partychange == 0,
            data = data) 
## nothing changes when you control for strength of partisanship or partisan social identity

## output this to a table 
models <- list()

models[["Behavior (1)"]] <- m2a 
models[["Behavior (2)"]] <- m2b 
models[["Face Mask (1)"]] <- m4a 
models[["Face Mask (2)"]] <- m4b 
models[["Safety (1)"]] <- m1a 
models[["Safety (2)"]] <- m1b 
models[["Policy (1)"]] <- m3a 
models[["Policy (2)"]] <- m3b 

## add the N 
rows <- tribble(~term, ~"m1a",~"m1b",~"m2a",~"m2b",~"m3a",~"m3b",~"m4a",~"m4b",
                'N',m1a$N,m1b$N,m2a$N,m2b$N,m3a$N,m3b$N,m4a$N,m4b$N)
attr(rows,'position') <- 30 

## clean up the GOF output (just N & R-Squared)
gm <- modelsummary::gof_map 
gm$omit <- TRUE 
gm$omit[5] <- FALSE 

msummary(models,
         add_rows = rows, 
         coef_map = c("apstable" = "Animus",
                      "rep" = "Republican",
                      "year" = "Wave", 
                      "pid_strength" = "Partisan Strength",
                      "sociden_pid" = "Partisan Social Identity",
                      "year:rep" = "Republican*Wave",
                      "year:apstable" = "Animosity*Wave",
                      "apstable:rep" = "Animosity*Republican", 
                      "year:pid_strength" = "Strength*Wave",
                      "year:sociden_pid" = "Soc. Identity*Wave",
                      "year:rep:pid_strength" = "Strength*Republican*Wave",
                      "year:rep:sociden_pid" = "Soc. Identity*Republican*Wave",
                      "year:apstable:rep" = "Animosity*Republican*Wave"), 
         gof_map = gm, 
         stars = T, 
         output =  "Tables/raw/Table_A5_2.html") 